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ABSTRACT 

Using the most recent measurements of the ultraviolet (UV) luminosity functions (LFs) 
and dust estimates of early galaxies, we derive updated dust-corrected star-formation 
rate functions (SFRFs) at z ~ 4 — 8, which we model to predict the evolution to higher 
redshifts, z > 8. We employ abundance matching techniques to calibrate a relation 
between galaxy star formation rate (SFR) and host halo mass Mh by mapping the 
shape of the observed SFRFs at 2 ~ 4-8 to that of the halo mass function. The resulting 
scaling law remains roughly constant over this redshift range. We apply the average 
SFR — Mh relation to reproduce the observed SFR functions at 4 < z < 8 and also 
derive the expected UV LFs at higher redshifts. At 2 ~ 9 and z ~ 10 these model LFs 
are in excellent agreement with current observed estimates. Our predicted number 
densities and UV LFs at z > 10 indicate that JWST will be able to detect galaxies 
out to 0 ~ 15 with an extensive treasury sized program. We also derive the redshift 
evolution of the star formation rate density and associated reionization history by 
galaxies. Models which integrate down to the current HUDF12/XDF detection limit 
{M-uv ~ -17.7 mag) result in a SFRD that declines as (1 + ; r )~ 10 ' 4±0 ' 3 at high redshift 
and fail to reproduce the observed CMB electron scattering optical depth, r ~ 0.066, 
to within ler. On the other hand, we find that the inclusion of galaxies with SFRs well 
below the current detection limit (Mjjy < -5.7 mag) leads to a fully reionized universe 
by z ~ 6.5 and an optical depth of r ~ 0.054, consistent with the recently derived 
Planck value at the ler level. 

Key words: cosmology: theory - galaxies: high-redshift - galaxies: mass function. 


1 INTRODUCTION 

The sensitive, near-infrared imaging capabilities of the Hub¬ 
ble Space Telescope (HST) have significantly expanded our 
understanding of galaxy evolution through cosmic time. Ef¬ 
forts to identify galaxies through photometric selections and 
follow-up spectroscopy have progressively extended the ob¬ 
servational frontier to higher and higher redshifts, resulting 


in robust galaxy samples in the range 2 ~ 4-8 (e.g. Madau 

et al.| 

1996. Steidel et al. 1999| |Bouwens et al. 2004[ 

2010 

2011a 

2014a 

Dickinson et al. 2004; Yan & Windhorst 

2004 

Bunker et al 

2010 Finkelstein et al. 2010, 2012a|b| 

2014 

McLure et al 

2010, 2013; Oesch et al. 

2010, 2012a| Schenker 

et al. 2012 2 

313 Bradley et al. 2014 

Schmidt et al. j2014). 


With the Wide Field Camera 3 (WFC3/IR) on board the 
HST, the current frontier for identifying high-redshift galax¬ 
ies now lies at z ~ 11, a mere half a billion years after the 
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Big Bang, and new surveys are now building up the sample 


sizes of 2 > 9 galaxies (e.g. Bouwens et al.||2011b ' 

311is et 

al. 2013; [Zheng et al. 2012,,Coe et al. 2013; McLure et al. 

2013; Oesch et al. 2013, 2014| Bouwens et al. 2014b; 

Oesch 

et al.|2012b McLeod et al. 2015 Ishigaki et al.|2015 

■ 


However, due to the small number statistics, our un¬ 
derstanding of galaxies beyond 2 ~ 8 is still very limited 
and obtaining accurate constraints on the evolution of the 
galaxy rest-frame ultraviolet luminosity function (hereafter 
referred to as the luminosity function, or UV LF) in the early 
universe remains challenging. Typically parameterized by a 
Schechter function with a power-law slope at faint luminosi¬ 
ties and an exponentially declining form at the bright end, 
the luminosity function provides a measure of the relative 
space density of galaxies over a wide range of luminosities 
at a particular redshift (Sc hechter|1976 ). Since the UV 1500 
A light probes recent star formation activity, the integral of 
the dust-corrected (intrinsic) luminosity function can then 
be used to derive the cosmic star formation rate density 
(SFRD) and estimate the galaxy population’s contribution 




















































2 N. Mashian et al. 


to cosmic reionization, a process that appears to have been 
completed by 2 ~ 6 (Finkelstein et al. 2012a Loeb & Furlan 


etto 


2013 Robertson et al. 2013 Bouwens et al. 2015a). 


Given the scarce number of 2 ~ 9-10 galaxy candidate de¬ 
tections, a reliable determination of the SFRD at these red- 
shifts remains difficult, with conclusions from separate anal¬ 


yses disagreeing on the SFRD evolutionary trends (Coe et 
al.|2013 Ellis et al.|2013 Oesch et al.|2013 20141. Further¬ 


more, with the identification of galaxies at these redshifts, 
we are approaching the limit of HST’s detection capabilities. 
Sources beyond 2 ~ 11 will continue to remain largely inac¬ 
cessible until the advent of the James Webb Space Telescope 
(JWST), whose scheduled launch in 2018 promises to yield 
revolutionary insights into early galaxy formation at 2 > 10 
( Gardner et al.|2006 ). 

The analysis and extrapolation of the UV luminosity 
function evolution towards higher (z « 8) redshifts has been 
the focus of several past studies, in the form of both hydro- 


dynamical simulations (e.g. Nagaminc et al. 


Finlator 


& Ozel 2011. Salvaterra, Ferrara. & Daval 2011 

Jaacks et 

al. 

2012) Dayal et al. 

2013 OShea et al. 2015) 

and semi- 

analytical models (e.g. |Trenti et al. 2010 

Munoz|2012 Tac- 

chclla, Trenti, & Carollo|2013; Cai et al. 

2014; Dayal et al. 


2014| Behroozi fe Silk||2015 l. Successfully reproducing the 
statistical properties of the observed Lyman-break galaxy 
(LBG) population at 2 = 6-9, simulations predict a large 
number of undetected, low-mass galaxies that may have sig¬ 
nificantly contributed to the reionization of the Universe 
at 2 > 6. Semi-empirical models which tie the evolution 
of galaxy luminosity to the dark matter properties of the 
host halos reach similar conclusions, predicting a steepening 
faint-end slope of the LF with redshift and a sharp decline 


in the cosmic star formation rate above 2 =8 (Trac, Ccn, & 
Mansheld|2015 Mason, Trenti, fe Treu|2015 l. 


In this paper, we take advantage of the recent progress 
in the observed galaxy statistics at redshifts 2 ~ 4-8 to 
evolve the luminosity function, along with its derived prop¬ 
erties, towards higher redshifts by establishing a link be¬ 
tween the galaxy star formation rate (SFR) and its host 
halo mass via abundance matching techniques (Vale & Os- 
|triker||2004| |2006| |Conroy fc Wechsler||2009[ ). Since galaxy 
formation is governed by the inflow of baryonic matter into 
the gravitational potential wells of dark matter halos where 
it cools and initiates star formation, the properties of galax¬ 
ies are invariably linked to the characteristics of their host 
halos. We thus seek to derive a scaling relation between dark 
matter halo mass, Mu and galaxy SFR at each redshift by 
assuming that there is a one-to-one, monotonic correspon¬ 
dence between these two properties. We anchor our model 
to the observed luminosity functions at 2 ~ 4-8 corrected 
for dust-extinction, mapping their shape to that of the halo 
mass function at the respective redshifts. Finding that the 
SFR-Mh scaling law remains roughly constant across this 
redshift range, we apply the average relation, SFR av (Mh), 
to reproduce the observed 2 ~ 9-10 luminosity functions and 
then further explore how the UV LF, and the corresponding 
SFRD, evolve at higher redshifts, 2 ~ 11-20. 

We describe the details of our approach and the theo¬ 
retical framework used to derive the shape and amplitude 
of the galaxy UV LF in §2. Our predictions for the evolv¬ 
ing LF, the resulting SFRD, and the contribution of these 
galaxy populations to cosmic reionization at high redshifts 


are presented in §3. We conclude in §4 with a summary of 
our findings and their implications for future surveys with 
JWST. We adopt a flat, ACDM cosmology with Q m = 0.3, 
Da = 0.7, = 0.045, Ho = 70 km s -1 Mpc -1 , i.e. h = 

0.7, ag = 0.82, and n s = 0.95, consistent with the most 


recent measurements from Planck (Planck Colla boration et 
ah|2015a|. All magnitudes in this work are in the AB scale 


( Oke fe Gunn|[l983 1. 


2 THE FORMALISM 

The goal of this paper is to derive an empirical prediction 
of the evolution of UV LFs at 2 > 8 tied to the evolution of 
dark matter haloes. We do so by first calibrating a SFR—Mh 
relation at 2 ~ 4 — 8, which we then evolve to 2 > 8. Our 
formalism is described in the following sections. 


2.1 The Observed Star-Formation Rate Functions 


One important limitation of the UV LFs in probing the 
galaxy build-up in the early universe is dust extinction, 
which significantly affects the UV luminosities of galaxies. 
Thanks to deep multi-wavelength photometry with HST, it 
has become clear that dust extinction evolves with redshift, 
becoming significantly less important at earlier cosmic times 


(Bouwenset^Jj|2009jjCastellano^t_alJ[2012jjFinkelstein_et 


al.|2012a |Dunlop et ai.|2013||Bouwens et al.|2014a||Rogers 


et al.||2014 1. Thus, the UV LF at high redshift evolves due 
to two effects: (i) evolution in dust extinction, and (ii) SFR 
build-up of the underlying galaxy population. Any model of 
the UV LFs therefore has to disentangle these two effects. 

Alternatively, we can directly make use of the SFR func¬ 
tions (SFRF), (f>(SFR , 2 ). These provide the number density 
of galaxies with a given SFR, and are thus corrected for dust 


extinction. Following Smit et al. (20121, we derive the ob¬ 


served SFR functions based on the most recent UV LFs at 


2 > 4 taken from Bouwens et al. (2015aI, which we correct 


for dust extinction using the observed LTV continuum slopes 
/3 and their relation with UV extinction ( Meurer et al.| 1 999j ). 
The UV continuum slopes are modeled as a function of lu¬ 
minosity using the most recent relations from Bouwens et 


. For more details on this approach see Smit et 


The dust correction effectively shifts the LF to higher 
luminosities and slightly lowers the volume densities due 
to the renormalization of the magnitude bins. The dust- 
corrected, intrinsic UV luminosities are then converted to 
SFRs using the following empirical relation (Kennicutt| 
19981), 


SFR 

M 0 yr" 


= 1.25 x 10" 


L 


UV,corr 


erg s 


^z" 


and the resulting SFR functions are represented with 
typical Schechter parameterization, 


(1) 

the 


<KSFR)dSFR 


SFR \“ 

sfitJ exp 


SFR \ d SFR 
SFR* ) SFR* ' *■ ’ 


The Schechter parameters for the 2 ~ 4-8 SFR functions 
are directly related to the Schechter function parameters 
of the UV LF and the slope of the Luv — fi relation (see 
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Table 1 . Schechter Parameters Determined for the SFR Functions 


<*> logl0 ^f? T 


3.8 

-2.79i0.07 

1.61i0.06 

-1.53i0.03 

4.9 

-3.25i0.10 

1.75i0.09 

-1.59i0.08 

5.9 

-3.45i0.16 

1.62i0.14 

-1.62i0.08 

6.8 

-3.67i0.23 

1.54i0.20 

-1.76i0.12 

7.9 

-3.79i0.31 

1.31i0.36 

-1.79i0.18 


sections 2.1 and 2.2 in Sm it et al.| |20l2j). Using the most 
recent UV LFs and UV continuum slope measurements, we 
thus update the previous SFR functions of |Smit et ah ( |2012[ ). 
The resulting Schechter function parameters are listed in 
Table [l] They are also shown in Figure [I] where we compare 
our results with the previous analysis of Smit et al. (20121. 
As can be seen, the combination of less evolution in the 
characteristic luminosity found in newer UV LF results, as 
well as the consistently redder UV continuum slopes found in 
updated measurements, result in significantly higher SFRFs 
compared to the previous analysis, in particular at higher 
redshift. 

We perform the same analysis also for z ~ 9—10 galaxy 
SFR functions. These are much more uncertain, given the 
very large uncertainties in the UV LFs at these redshifts. 
Also, note that based on the evolution of the UV continuum 


slope distribution at lower redshift (Bouwens et al. 2012a 


2014a[ |Dunlop et al.|2013[ |Finkelstein et al.||2012a| |Wilkins| 


et al.|20~ i, the dust correction is expected to be negligible 

at redshifts z > 8; the z ~ 9 and 10 SFR function parameters 
are thus directly related to the uncorrected UV LFs at those 
redshifts. 


2.2 Method to derive the average SFR — Mh 
relation 


Adopting the approach presented in Vale & Ostriker (2004), 


we use the observed, dust-corrected z ~ 4-8 luminosity func¬ 
tions, i.e. SFR functions, to derive an empirical relation be¬ 
tween the galaxy SFR and its host halo mass. In this abun¬ 
dance matching method, the SFR-Mj, relation is calculated 
by setting the SFR of a galaxy hosted in a halo of mass Mh 
to be such that the number of galaxies with a star formation 
rate greater than SFR equals the number of halos with mass 
greater than Mh at a given epoch: 


r oo 

/ 

Jm h 


r oo 

n(M' h , z)dM' h = 0(SFR’, z)dSFR\ (3) 

J SFR 


where the </>(SFR, z) are the observed SFR function derived 
in the previous section and n(Mh,z) is the Sheth-Tormen 
halo mass function (Sheth & Tormen 19991, 


n(Mh, z)dM = A 



/FprndV- _„ 2/2 

\ tv M dM 


(4) 


with v — y/aS c /[D(z)a(M )], a = 0.707, A = 0.322 and q 
= 0.3; <j(M) is the variance on the mass scale M (assum¬ 
ing the linear theory density power spectrum) while D(z) is 
the growth favor and S c is the linear threshold for spherical 
collapse, which in a flat universe is S c = 1.686. 



log SFR (M sun /yr) 

Figure 1. The evolution of the star-formation rate function, 
‘f’SFR) across redshift z = 4 — 8. The points correspond to the 
dust-corrected stepwise UV LFs from |Bouwens et al.] ( |2015a[ ), 
which were corrected for dust extinction based on the UV con¬ 
tinuum slope distributions measured in [ Bouwens et al.| ( [2014a ). 
The solid lines show the corresponding Schechter parameters for 
the SFRFs. These can be compared to the dashed lines from|Smit| 
|et al.| ( [20l2t , which were derived using the same procedure, but 
based on now-outdated UV LFs and /3 distributions (covering 
only z = 4 — 7). The combination of less evolution in the charac¬ 
teristic luminosity found in newer UV LF results as well as the 
consistently redder UV continuum slopes, result in significantly 
higher SFRFs compared to the previous analysis, in particular at 
higher redshift. 


The abundance matching technique, as prescribed in 
equation (3), presupposes that the SFR is a monotonic func¬ 
tion of the halo mass. This is supported by the observed 
trend that the clustering strength of galaxies at high red¬ 
shift increases with their UV luminosity, similar to that of 


the clustering strength of halos increasing with mass ( Lee et 
|al.|2009[ ). Hence, a scaling relation that assumes the SFR of a 
galaxy is a monotonically increasing function of the mass of 
its host halo is a reasonable starting point. Furthermore, our 
model assumes a one-to-one correspondence between galax¬ 
ies and host halos. This neglects the substructure expected 
to exist within a halo, comprised mainly of halos that formed 
at earlier epochs and merged to become subhaloes of more 
massive hosts (Bullock, Wechsler, & Somerville 2002 Wech- 


sler et al.|2000 1. To include the contribution of this subhalo 

population and the galaxies they are potentially hosting, we 
modified the left-hand side of equation (3) to integrate over 
the total mass function, i.e. the sum of the regular Sheth- 
Tormen halo mass function and the unevolved subhalo mass 


function taken from Vale & Ostriker (2004). We found that 


because values for the subhalo mass function are generally 
lower than those of the halo mass function, the inclusion of 
this additional substructure has a negligible effect on the 
resulting scaling law; the subhaloes’ contribution to the ex¬ 
pression is most important at the low-mass end, and even 
there, its inclusion changes the results by < 10%. Therefore, 
for the sake of simplicity, we neglect multiple halo occupa¬ 
tion and assume that each dark-matter halo hosts a single 
galaxy. 
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Figure 2. Top panel: the relation between galaxy SFR and dark matter halo mass, SFR(Mf l , z), derived at z ~ 4-8 via abundance 
matching of the dust-corrected SFR functions to the halo mass function. For comparison, we show the SFR — M relation at redshift z = 
8 obtained by abundance matching the observed (uncorrected) UV LF function to the halo mass function at z = 8 (dotted black line). 
The solid black curve represents the mean SFR-A//, scaling law, SFR av (A'/), ), obtained by averaging the (dust-corrected) relations over 
the redshift range at a fixed halo mass. The uncertainty in our model calibration of SFR av (M (shaded region), was derived by varying 
the z ~ 4-8 SFR function parameters within their respective la confidence regions. Bottom panel: ratio of the mean scaling law to the 
SFR-AT), relations derived at each redshift. The vertical axis represents the factor by which the average relation over/underestimates the 
z ~ 4-8 scaling laws over the halo-mass range of interest. 


2.3 Modeling the SFR Functions 


Once the average relation between the SFR and the halo 
mass are known, we can model the SFR functions assuming 
a realistic dispersion around the average relation. We model 
the probability density for a halo of mass Mu to host a 
galaxy with a star formation rate SFR to obey a log-normal 
distribution (Giavalisco &: Dickinson|2001 Yang, Mo, & van 


den Bosch|2003 Munoz fe Loeb||2011l, 


P(SFR\M h )dSFR = 

1 ( log? 0 (SFR/SFR av (M h ))\ J , OT ,„ /c , 

^ eX H-2^-JrflogtoSFR (5) 

where SFR av (Mh) is the mean SFR-M^ relation derived by 
averaging the SFR(Mh, z) relations obtained for 4 < z < 
8 over the specified range of redshifts at fixed halo mass. 
The variance, a 2 , represents the scatter in SFR at a fixed 
Mh that arises from the stochastic nature of star forma¬ 
tion activity. Halos of similar masses can have a range of 
large-scale environments, merger histories, and central con¬ 
centrations (Lee et al. '2009). These conditions, along with 
interactions with nearby systems can result in different rates 
of gas accretion and star formation in galaxies. To account 
for this variance in the SFR-M/j relation, we adopt a con¬ 
stant intrinsic scatter of 0.5 dex. This choice is motivated by 


the 0.5 dex scatter in stellar mass at constant SFR found by 


Gonzalez et al. 

(2011), a result further confirmed by Wyithe, 

Loeb, & Oesch 

(20141 as providing the best statistical fit to 


observations of the SFR-TVL, relation. 


The SFR function, ij>{SFR), can then be obtained by 
integrating over the probability-weighted number densities 
of all halos that can achieve the star formation rate SFR 
within the allowed scatter: 


<j>( SFR, z) = 


/ 


dM h n{M h ,z)P{SFR\M h ) 


1 1 
aV^SFR 


log? n (SFR/SFR. a „(M h )) 

dM h n(M h ,z)e 2 » 2 


( 6 ) 


Note that in the limit of the variance a 2 —> 0, the log-normal 
probability distribution becomes a delta function and the 
equation is reduced to the simpler form, 


0(SFR, z) = n(M h , z) 


dSFR a „(M h ) 

dM h 



(7) 
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Figure 3. Comparison of the predicted SFRFs at z ~ 4 — 8 assuming an intrinsic scatter a = 0.5 dex (equation (6); solid colored curves) 
with those obtained by setting the variance equal to zero (equation (7); dashed curves). Our model SFRFs are overplotted with the 
data points at each respective redshift, along with the Schechter functions found to best approximate these measurements (gray curves; 
parameters recorded in Table 1). The data points represent the intrinsic UV luminosities, corrected for dust extinction and converted to 
SFRs using equation (1); the SFR functions are thus the intrinsic, dust-corrected LFs at these redshifts. The shaded areas in each panel 
correspond to the lcr confidence regions for our model LFs (when intrinsic scatter is accounted for). 


3 RESULTS 


3.1 The Constant SFR—Mh Relation 

The curves in the top panel of Figure [2] depict the SFR-M^ 
relations inferred from the observed LF functions at redshifts 
4 < z < 8 via the abundance matching technique discussed 
in §2.1. The scaling law remains roughly constant within 0.2 
dex across this redshift range and is fairly well fit by a double 
power law with a turnover at a characteristic halo mass of 
as 2x 10 11 Mg. Averaging the relations over z ~ 4-8 at a 
fixed halo mass, we obtain a mean scaling law, SFR av {Mh) 
(solid black curve), characterized by a power-law slope /3 ~ 
0.9 at the high-mass end, i.e. Mh > M£, and /3 ~ 1.5 at the 
faint end where Mh falls below the characteristic mass. 

These results are consistent with the models presented 


in Trenti et al. (2010) where the L — Mu relation is calibrated 
at z ~ 6 and ai~ M^ relation is derived with /3 ~ 1.3-1.6 
at the low-mass end. This steepening of the SFR-Mj, relation 
towards lower masses may be due to the increased impor¬ 
tance of feedback mechanisms in the star formation activity, 


such as supernovae and stellar winds (Giavalisco & Dickin¬ 
son pOOl). Meanwhile, the overall flattening of SFR(Mh) at 
these redshifts compared to the local universe, where SFR 
oc is mainly a consequence of the steepening faint-end 
slope of the observed LF at these redshifts, where a < -1.5 
( Cooray fe Milosavljevic] 2005). 

The uncertainty in our model 
SFR av {M h ), represented by the shaded 
ure [2j is derived by varying the z ~ 4-8 SFR function 
parameters within their respective la confidence regions 
and rederiving the average SFR-M^ relation. As can be 
seen, the dashed curves defining the SFR-M^ relation 
at each redshift fall well within the boundaries of the 


calibration of 
region in Fig- 


uncertainty in the average scaling law, with the exception of 
the z ~ 4 relation at the high-mass end. The bottom panel 
of Figure [2] illustrates the accuracy with which the mean 
scaling law approximates SFR as a function of halo mass 
at different redshifts. Overall, the SFR av {Mh ) relation 
provides a reasonable representation of the SFR(Mh,z) 
relations at all redshifts, with departures of less than 
0.2 dex within the halo mass range of interest. With the 
exception of the z ~ 4 relation, the SFR(Mh, z) are all 
very similar (within 0.5 dex) at Mh > 10 12 Mg. At the 
low-mass end, the z ~ 5-6 and z ~ 7-8 scaling laws are 
over- and under-estimated respectively by less than a 
factor of two, while the SFR at z ~ 4 is overestimated 
over the full halo-mass range by approximately the same 
factor. We therefore expect the mean relation to yield 
accurate estimates of the SFR function at these redshifts 
and, assuming that the SFR-Aifi, relationship continues to 
remain roughly unchanged for z > 8, at higher redshifts as 
well. 


3.2 The Modeled SFR Functions at z = 4 — 8 

Applying SFR av (M h ), we derive the expected SFR func¬ 
tions at redshifts z > 4 assuming only evolution of the un¬ 
derlying dark-matter halo mass function. The results of our 
LF model are shown in Figure [3] both for the case where 
an intrinsic scatter of a ~ 0.5 dex is assumed (solid col¬ 
ored curves) and in the limiting case where the variance is 
set to zero (dashed curves). Plotted alongside these results 
are the dust-corrected data points and the Schechter LFs 
with the observed best-fit parameters reported in Table 1 
(gray curves). We find that including an intrinsic scatter in¬ 
creases the number density and yields SFR functions that 
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Figure 4. Comparison of the predicted LFs assuming an intrinsic scatter a = 0.5 dex (equation 6; solid colored curves) with those 
obtained by setting the variance equal to zero (equation 7; dashed curves). Since at £ ~ 9 and 10, the dust extinction is expected to be 
negligible, the curves represent the observed, uncorrected LFs. We include determinations of the z ~ 9 and z ~ 10 UV LFs from [Oesch| 
|et al.| ( |2013[ ) , |Oesch et al.| ( |2014| ) , and |Bouwens et al.| < |2015b[ ) . Arrows signify upper bounds on the number density at a given absolute 
UV magnitude, M u y . The shaded area in each panel represents the lcr confidence region for the model LF when intrinsic scatter is 
accounted for. Our model UV LFs are in excellent agreement with the current observed estimates at these redshifts (albeit the latter are 
still based on small samples). 


are more consistent with observation, particularly at the lu¬ 
minous end where the increase in <f>(SFR) due to scatter is 
most apparent. This boost in number density with the in¬ 
troduction of scatter can be attributed to the shape of the 
halo mass function. Even though SFR(Mh) is log-normally 
distributed around the mean, SFR av (Mh), and the SFR 
scatter can thus go in either direction, the net change in the 
LF will always be dominated by low-mass halos entering into 
the galaxy sample by scattering into a SFR higher than its 
mean value, leading to a larger estimate of rf)(SFR). 

With the exception of 2 ~ 4, our predicted SFR func¬ 
tions accurately reproduce the observed ones and provide 
very good fits to the data points at the luminous end when 
intrinsic scatter is accounted for. At z ~ 4 our model pre¬ 
dicts a SFRF which lies significantly above the observed 
data points by factors of 2-10 x. This is likely due to a com¬ 


bination of two effects: ( i ) As also argued by Smit et al. 
|2012]), the bright end of the observed z ~ 4 SFRF is likely 
underestimated due to heavily obscured galaxies, which are 
missed in UV-selected Lyman break galaxy samples; and (ii) 
quiescent galaxies start to appear in the Universe at z ~ 4 
in significant numbers (e.g. Straatman et al. 2014). Since 
our model does not include any prescription to quench star 
formation in galaxies, this is likely another reason why the 
model appears to be inadequate at z ~ 4. We conclude that 
our model should only be used at z > 4, where it provides 
a very good representation of the observed SFRFs. 


3.3 A Prediction for the UV LF at z = 9 — 10 

Building on the successful calibration of the SFR-M/j rela¬ 
tion at z ~ 5 — 8, we use our model to predict the galaxy 


UV luminosity function at z ~ 9 — 10. This represents the 
most distant galaxies in reach with HST, and recent observa¬ 
tions have started to provide the first LF estimates at these 
redshifts. However, the sample sizes are still very 
suiting in highly uncertain LF parameters (see e.g. 
aLp0 13l. 


small re- 
Oesch et 


We derive the model UV LFs at z = 9 and 10 by convert¬ 
ing the SFR -Mh relation calibrated in the previous section 
to a Luv ~ Mh relation using the standard Kennicutt (19981 
SFR-L(jy relations. As can be seen in Figure |4| our model 
provides an excellent fit to the current, albeit scarce, data. 
In particular, our model explains the somewhat puzzling ab¬ 
sence of z ~ 10 galaxy candidates between Muv = —18.5 
and —20.5 in current datasets. Incoming observations of the 
Hubble Frontier Field initiative will allow us to test this 
model further at z ~ 9 — 10 in the future. 


3.4 Extrapolation to z ~ 20 and Predictions for 
JWST 

Motivated by the good agreement obtained by applying an 
unevolving SFR -Mh relation to reproduce the observed LFs 
at 5 < z < 10, we extend our approach to even higher red¬ 
shifts for JWST predictions. The predicted LFs for redshifts 
z > 10, again assuming a constant intrinsic scatter of a = 
0.5 dex, are shown in Figure[5] These model LFs, which are 
well-fitted by Schechter functions across the redshift range 
4 < z < 20, significantly evolve with redshift. The drop in 
number density from z ~ 4 to 20 is reflected in the evolu¬ 
tion of the characteristic number density, <j >*, which decrease 
as d\og(f>*/dz ~ -0.3. We also find a gradual steepening of 
the faint-end slope ( da/dz ~ —0.08) and a shifting of the 
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Figure 5. Predicted LFs for 2 : > 10 based on our empirical model. 
The shaded region indicates the observable magnitude and vol¬ 
ume within reach of a 200 arcmin 2 JWST treasury sized ultra¬ 
deep survey (mn m — 31.5 mag). This indicates that the limiting 
redshift range observable with JWST may be 2 : ~ 15 for unlensed 
galaxies. 


characteristic luminosity towards fainter values for increas¬ 
ing redshift (dM^y/dz ~ 0.4). These trends in Schechter 
parameters are quite consistent with results from previous 
studies and extrapolations of observed z < 10 LFs, except 
that our model M* evolves more significantly than found in 
observ ations (e.g., |Trenti et al. 2010. Bo uwens et al.pO lla} 


2012a McLure et al.|2013 Finkelstein et al.|2014l. 

The shaded region in Figure |5] corresponds to the co¬ 
moving volume and magnitude range accessible with a 200 
arcmin 2 treasury sized ultra-deep JWST survey (m/ lm ~ 
31.5 mag). According to our model predictions, the limiting 
redshift for galaxy observations with a very deep JWST sur¬ 
vey is z ~ 15 (see also Windhorst et al.[2006). At higher red- 


shifts, where the surface density drops below ~10 -6 Mpc~ 3 , 
wider, deeper surveys would be required to observe the lu¬ 
minosity function of these rare, faint objects. 

Figure [6] shows the model predictions for the redshift 
evolution of the SFR density, obtained by integrating the 
SFR functions down to different SFR limits, ranging from 
SFRmin ~ 10 -5 Mo/yr to 0.7 Mg/yr. A minimum SFR of 
0.7 Mg/yr corresponds to Muv ~ -17.7 mag, the magni¬ 
tude of the faintest object observed in the HUDF12/XDF 
data; integrating down to this limit thus facilities compar¬ 
ison with the most recent, dust-corrected measurements of 
the SFRD, which have been plotted alongside the predicted 
curves in Figure [6] While our model appears to overpredict 
the SFRD measured for z < 6 by ~ 0.2-0.3 dex, the esti¬ 
mates of the p* at higher redshifts agree quite well with ob¬ 
servations. Furthermore, the evolution of the cosmic SFRD 
in our model is consistent with previous published results: 
the dust-corrected SFRD values evolve as (1 + z j~ 4 - 3±01 
at 3 < 2 < 8 before rapidly declining at higher redshifts 
where p* oc (1 + z )- 10 - 4±0 - 3 for 8 < z < 10. In addition 
to accurately reproducing the order of magnitude drop in 
SFRD from z < 8 to 2 < 10 that has been observationally 
inferred in several separate analyses ( |Oesch et al.]2013j), our 
model finds that the SFRD continues to steeply decline to- 



sity (SFRD) p* derived by integrating the model SFR, functions 
down to different star-formation limits, ranging from SFR m i n ~ 
10 -5 to 0.7 Mg/yr. The shaded blue and red regions repre¬ 
sents the la uncertainty in the SFRD obtained when integrating 
down to M jjy ~ -17.7 mag (i.e., SFR m i n ~ 0.7 Mg/yr) and 
Muv ~ -5-7 mag (i.e., SFR m i n ~ 10 -5 Mg/yr), respectively. 
Our model results (blue curve) can thus be compared with the 
(dust-corrected) measurements of the cosmic SFRD (blue circles) 
derived from this dataset. We find that the best-fit evolution of 
the SFRD at 8 < 2 < 10 is p* oc (l+z) _10 - 4±0 - 3 , significantly 
steeper than the lower redshift trends which fall as (l+z) —4 3 . 


wards higher redshift, predicting a cosmic SFRD of logp* ~ 
-7.0±0.3 Mg yr _1 Mpc” 3 at z ~ 16. 

However, as will be shown and discussed in the follow¬ 
ing section, a star-formation rate density that declines as 
(1 + z) -10 ' 4 at high redshifts, as derived assuming a min¬ 
imum SFR of 0.7 Mg/yr, fails to reproduce the observed 
Planck optical depth. Furthermore, galaxies are expected to 
exist beyond the current, observed magnitude limit. The¬ 
oretical and numerical investigations indicate that a halo 
at z < 10 irradiated by a UV field comparable to the one 
required for reionization needs a minimum mass of Mu ~ 
6x10' Mg in order to cool and form stars 

. We therefore explore the 
integrating the SFR den¬ 
sity down to minimum SFRs as low as io _ l 5 0±0 - 8 ) Mg yr _1 , 
the star formation rate corresponding to the minimum halo 
mass based on our average SFR-M^ relation. This leads to a 
more moderate decline of the cosmic SFRD towards higher 
redshift. In this case, the best-fit evolution for 8 < z < 10 is 
p* oc (1 + z) _ ' ■ 3±0 ' 5 , with an estimated SFRD of ~ 4xl0 -5 
Mg yr -1 Mpc -3 at z ~ 16, a mere ~ 250 million years after 
the Big Bang. 


& Loeb 1996 Tegmark et al.|1997 


implications of this prediction by 


Haiman, Thoul, 


3.5 Contribution of Galaxies to Cosmic 
Reionization 

If star-forming galaxies supply the bulk of the photons that 
drive the cosmic reionization process, then our redshift- 
evolving SFRD p*(z) can be used to determine the reion¬ 
ization history of the universe. The ionized hydrogen frac¬ 
tion Q(z) can be expressed as a time-dependent differential 
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Figure 7. Left panel: The reionization history Q(z) calculated by solving the differential equation given by equation (8) assuming an 
ionizing photon escape fraction f esc = 0.2, an IGM clumping factor of C = 3, and an average number of ionizing photons per unit mass 
of stars Nry = 7.65xl0 60 Mq 1 . The blue and red curves denote the reionization histories derived by using the cosmic SFRD p* integrated 
down to SFR rn i n = 0.7 and 10 — ( 5 - 0 ± 0 - 8 ) M©/yr respectively (i.e. M jjy < -17.7 and < -5.7 mag; galaxies hosted by the minimum halo 
mass to cool and form stars). For these chosen parameters, 68% confidence limits on the redshift of half-reionization are 6.2 < z < 7.3 
when considering only the currently observed galaxy population, and 6.5 < z < 8.1 with the inclusion of galaxies that expected to exist 
beyond this magnitude limit. Right panel: The corresponding Thomson electron scattering optical depths r integrated over redshift from 
the present day, along with the Planck constraint r = 0.066d=0.016 (gray area). The shaded regions correspond to the 68% confidence 
intervals, computed from the uncertainties in the cosmic SFRD p* shown in Figure [6] 


equation (jLoeb & F urlanetto|2013 ) 


Q = 




Q 


(n H ) 




( 8 ) 


where (uh) is the comoving number density of hydrogen 
atoms and riion is the comoving production rate of ioniz¬ 
ing photons, hion = fescN 1 p,(z) where IV 7 is the number of 
ionizing photons per unit mass of stars and f esc is the aver¬ 
age escape fraction from galaxies. The recombination time 
in the IGM is 


tree — \CcXb{T^)1 1 -(- Y p /A:Xp) (jlH ) (1 + z) ] 


(9) 


where afl(T) is the case B recombination coefficient for hy¬ 
drogen (we assume an IGM temperature of 10,000 K cor¬ 
responding to olb ~ 2.79xl0~' 9 Mpc 3 yr _1 ), X p = 0.75 
and Y p = 1-X P are the primordial hydrogen and helium 
mass-fractions respectively, and C is the clumping factor 
that accounts for the effects of IGM inhomogeneity. Esti¬ 
mates of the contribution of star-forming galaxies to reion¬ 
ization therefore depend on assumptions about the stellar 
IMF, metallicity, escape fraction, and clumping factor. 

In this paper, we take JV 7 equal to 7.65 x 10 6 ° M7 1 , cor¬ 


responding to a low-metallicity Chabrier IMF (Behroozi & 
|Silk|2015| ). The escape fraction is more difficult to constrain, 
especially at higher redshifts where the correction for inter¬ 
vening IGM absorption systems is high; while inferences at 
z < 5 suggest fesc < 10%, simulations suggest that f esc can 


be very large (Wise & Cen 2009 Wise et al. 2014 Hayes 


et al. 2011) . For the sake of simplicity, we follow the ap¬ 


proach in the previous literature and adopt a constant fidu¬ 
cial value of fesc = 0.2 for all redshifts. The determination of 


the clumping factor is also uncertain, with estimates rang- 


ing from C = 2 to C = 4 at high redshifts ( 

Pawlik, Schaye, 

& van Scherpenzeel 2009. Finlator et al. 2012) Shull et al. 

2012 

1; we therefore assume C = 3 in our fiducial model (see 

also 

Oesch et al. 2009; Bouwens et al. 2012b; Finkelstein 

et a 

. 2014; Behroozi & Silk 2015. Bouwens et al. |2015c| 

Robertson et al. 2015 

I- 


For a given reionization history Q(z), the electron scat¬ 
tering optical depth r can also be calculated as a function 
of redshift, 


t(z) = / e(nrf}crT(l + Y p /4:X p )Q(z')H(z')(l + z') 2 dz' 

Jo 

( 10 ) 

where c is the speed of light, ot is the Thomson cross section, 
and H(z) is the Hubble parameter. This optical depth can 
be inferred from observations of the cosmic microwave back¬ 
ground, and for a while, the accepted value from the Wilkin¬ 
son Microwave Anisotropy Probe (WMAP) 9-year dataset 
was t = 0.088± 0.014, which, in the simplest model, corre¬ 
sponds to instantaneous reionization at z re ion — 10.5±1.1. 
In early 2015, a significantly lower value of r = 0.066±0.016 


was reported (Planck Collaboration et al. 2015aI, consis¬ 


tent with instantaneous reionization occurring at Zreion — 
8-8l 3 '^. More recently, the value has again shifted up to r 
= 0.078±0.019, consistent with instantaneous reionization 


occurring at Zreion — 9.91 3 '® (Planck Collaboration et al. 
|2015b[ ). 

Figure [7] shows the expected reionization history of the 
universe (left panel) and the corresponding optical depth 
(right panel) computed using the cosmic SFRD p* we de¬ 
rived from our model UV LFs. We find that the currently 
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observed galaxy population at magnitudes brighter than 
Muv < -17.7 (SFR > 0.7 M 0 /yr, blue curve) fully reionize 
the universe by redshift z ~ 6 and that the model corre¬ 
sponding to this magnitude limit predicts a Thomson opti¬ 
cal depth of r ~ 0.048. Integrating further down to Muv < 
-5.7 (SFR > 10 -5 M 0 /yr, red curve) produces a model 
that reionizes the universe slightly earlier, at z ~ 6.5, with 
a predicted optical depth of r ~ 0.054. The shaded regions 
corresponding to each of these curves represent the 68% con¬ 
fidence limits arising solely from the uncertainties in the 
normalization of the cosmic SFRD. 

Our results for the reionization history are consistent 
with measurements from high-redsliift quasar and gamma- 
ray burst spectra, as well as from Lyman alpha emission 
in high-z galaxies which collectively indicate that reioniza- 


& Faucher-Giguere 

2012| Chornock et al. 2013| Treu et al. 

2013 

Schenker et a 

.(2014. McGreer, Mesinger, L D’Odorico 


20151. Furthermore, the model which includes the popula¬ 
tion of low-mass halos with SFRs as faint as 10~ 5 M 0 /yr 
produces an optical depth that is consistent with the previ¬ 
ous inferred value by Planck, r = 0.066 ( |Planck Collabora7| 
tion et al.|2015a l, at the lcr level. These findings, in line with 
previous analyses, indicate that a significant population of 
low-mass star-forming galaxies is necessary for cosmic reion¬ 
ization and strengthen the conclusion that the bulk of pho¬ 
tons responsible for reionizing the early universe emerged 
from ultra-faint galaxies. 

In order for our model to reproduce the most recent 
value of r = 0.078 inferred from CMB observations by 
Planck, the average SFR — Mh relation used would have 
to allow for larger star formation rates at the low-mass 
end. One simple way to satisfy such conditions is to per¬ 
mit SFR av (Mh) to vary according to the average relation 
derived via abundance matching for Mh > 10 11 M 0 , but 
decline as a power law for Mh < 10 11 M 0 with a slope of 
P ~ 1.3, shallower than the slope derived empirically in §3.1 
at the low mass end. In such a model, the corresponding 
optical depth would be r ~ 0.076, similar to the latest re¬ 
duced measurement of r = 0.078. While the high-redshift 
LF functions derived assuming such a form for the average 
SFR— Mh relation are in tension with current observations, 
future measurements taken with JWST will shed more light 
on how the LF functions evolve at these low SFRs. 

We also note that the results for Q(z) and r(z) depend 
on the choices for the escape fraction / esc , the number of 
ionizing photons per unit stellar mass N-y, and the clump¬ 
ing factor C; if one assumes a smaller clumping factor, or 
a higher value for f esc or N 7 , the evolution of Q(z) will be 
shifted towards higher redshifts and an optical depth closer 
to the one measured by Planck Collaboration et al. ( |2015b I 
will be obtained. Furthermore, the ionizing photons of Pop¬ 
ulation III stars are expected to additionally increase the 
optical depth at high redshift (Sobral et al. 2015). 


4 SUMMARY 

In this paper, we use the most recent z ~ 4-8 UV LFs and 
UV continuum slope measurements to derive an empirical 
prediction of the evolution of UV LFs at z > 8. Assuming a 
monotonic, one-to-one correspondence between the observed 


galaxy SFRs and the host halo masses, we map the shape 
of the observed z ~ 4-8 UV LFs to that of the halo mass 
function at the respective redshifts. We find that the result¬ 
ing SFR — Mh scaling law remains roughly constant over 
this redshift range and is fairly well fit by a double power 
law, SFR oc Mh, with /3 ~ 0.9 at the high-mass end, i.e. 
Mh > 2xlO rl M 0 , and /3 ~ 1.5 at the low-mass end. We 
note that the unevolving nature of this relation with redshift 
is an empirical result, an explanation for which lies beyond 
the scope of this paper. Future work based on semi-analytic 
models and numerical simulations, as well as anticipated ob¬ 
servations with Atacama Large Millimeter Array (ALMA), 
are expected to shed further light on this matter. 

Applying this average SFR — Mh relation with an in¬ 
trinsic scatter of a ~ 0.5, we accurately reproduce the ob¬ 
served SFR functions at 5 < z < 10 (Figure 3 and 4) and 
extend our approach to predict the evolution of the UV LF 
at redshifts z ~ 11-20. We find an evolving characteristic 
number density <j>* which decreases as dlog</>*/dz ~ -0.3, a 
gradually steepening faint-end slope, da/dz ~ —0.08, and 
a shifting of the characteristic luminosity towards fainter 
values, dMuv/dz ~ 0.4. Given the comoving volume and 
magnitude range of an ultra-deep JWST survey, our model 
predicts that observations of the LF up to z ~ 15 are within 
reach in the absence of gravitational lensing, while deeper 
and wider surveys would be necessary to observe higher red¬ 
shift objects. 

We also derive the evolution of the SFR density by inte¬ 
grating the SFRFs down to various SFR limits, ranging from 
0.7 M 0 /yr (corresponding to the faintest object observed 
in HUDF12/XDF), to 10 -5 M 0 /yr (corresponding to the 
minimum halo mass necessary to cool and form stars). The 
inclusion of galaxies with SFRs well below the current de¬ 
tection limit results in a more moderately declining cosmic 
SFRD with redshift and leads to a fully reionized universe 
by z ~ 6.5. Furthermore, the corresponding predicted op¬ 
tical depth in this model, r ~ 0.054, is consistent with the 
reduced value of r = 0.066±0.016 inferred from CMB obser¬ 
vations by Planck at the lcr level. These results strengthen 
the claim that a significant portion of the reionizing pho¬ 
tons in the early universe were emitted by a population of 
low-mass, star-forming galaxies that have thus far evaded 
detection. 
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